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Abstract. 

Lagrangian Coherent Structures (LCS) act as the organizers of transport in fluid 
flows and are crucial to understand their stirring and mixing properties. In the case 
of oceanic flows, LCS are known to drive biological dynamics, from plankton to top 
predators, which stresses the importance of their characterization in realistic flows. 
Lyapunov exponents are useful tools to compute LCSs. In this paper we have used the 
Finite-Size Lyapunov Exponent (FSLE) to identify LCSs in two different types of three- 
dimensional turbulent velocity flelds. First, in a canonical turbulent flow (channel flow 
between two parallel plates) the LCSs have a complex three-dimensional shape and 
are advected by the flow. Second, in an oceanographic setting (a regional simulation 
of the Benguela area) the LCSs also show a complex pattern on the horizontal but the 
small vertical velocities typical of oceanic flows result in a curtain-like shape. 
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1. Introduction 

The use of stretching quantifiers such as the Lyapunov exponents, which measure 
the relative separation of transported particles [H EJ [3l H], has largely improved the 
Lagrangian analysis of fiuid fiows. On the one hand Lyapunov methods provide 
information on time scales for dispersion processes, with its relevance for mixing and 
stirring of fiuids jU El El EJ El [7j. On the other, they are useful to detect the so- 
called Lagrangian Coherent Structures (LCS). LCSs ^ are templates for particle 
advection in complex fiows, separating regions with difiFerent dynamical behavior and 
signalling the existence of barriers and avenues to transport, fronts or eddy boundaries 

The relationship of LCSs with ridges (local maxima) of Lyapunov fields has been 
soundly established for the case of finite-time Lyapunov exponents (FTLEs) [TH [15], 
although it should be mentioned that techniques more precise than Lyapunov methods 
are available |TT1 [16]. In our work, we use instead finite-size Lyapunov exponents 
(FSLEs), which quantify the separation rate of fiuid particles between two given distance 
thresholds [1^ 2j. They turn out to be convenient for the case of bounded fiows in which 
characteristic spatial scales are more direct to identify than temporal ones. Following 
many previous works ^IQl [6l [171 [121 HH] we assume that the mathematical results existing 
for FTLEs are valid for FSLEs to a good approximation. In particular we assume that 
LCSs can be computed as ridges of FSLEs, and that they are transported by the fiow as 
material surfaces/lines, with no fiux of particles through them. Observations presented 
here are consistent with those assumptions. 

Despite its relevance in real fiows, the full three-dimensional (3d) structure of 
LCSs is still an open subject. In 3d fiows, LCS were explored in atmospheric contexts 
121], and in a turbulent channel fiow at Rcr = 180 in l22j. A kinematic ABC fiow 
was studied in [23]. In the ocean, where it is widely recognized that filamental structures, 
eddies, and in general oceanic meso- and submeso-scale structures have a great infiuence 
on marine ecosystems [24l [25l [26l [27], the identification of LCSs and the study of their 
role in the transport of biogeochemical tracers has primarily been restricted to the two- 
dimensional (2d) surface layers. There are two concurrent reasons for this: a) because 
of stratification and rotation, vertical motions in the ocean are usually very small when 
compared to horizontal displacements; b) synoptic measurements (e.g. from satellites) of 
relevant quantities are restricted to the surface. A few previous results for Lagrangian 
eddies in 3d were obtained in Refs. pHl [29], by applying the methodology of lobe 
dynamics and the turnstile mechanism. Also, Refs. [30l [3l] used 3d FSLE fields to 
identify LCS in oceanic fiows. In particular, a mesoscale eddy in the Southern Atlantic 
was studied in and it was shown that oceanic LCSs presented a vertical curtain-like 
shape, i.e. they look mostly like vertical sheets, and that material transport into and 
out of the mesoscale eddy occurred through filamentary deformation of such structures. 

In this paper, we use 3d fields of FSLE to identify LCSs in a turbulent channel fiow 
and in an oceanic fiow. Observations of the similarities and differences between the two 
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systems, both in their computation and their physical meaning, helps to appreciate the 
power and scope of this Lagrangian technique in the analysis of fluid flows. In Section 
[2] we describe the methodology used to identify LCSs in 3d turbulent flows. Sections [3] 
and [4] are devoted to the turbulent channel flow and the oceanic flow, respectively, and 
Section [5] presents conclusions and ideas for future work. 



2. Methods 



2.1. Finite- Size Lyapunov Exponents 

In order to study non-asymptotic dispersion processes such as stretching at flnite scales 
and bounded domains, the flnite size Lyapunov Exponent was introduced [UEIIS]. It is 
deflned as: 

A = -log^, (1) 
T do 

where r is the time it takes for the separation between two particles, initially rfo, to reach 
a value df. In addition to the dependence on the values of do and rf/, the FSLE depends 
also on the initial position of the particles and on the time of deployment. Locations 
(i.e. initial positions) leading to high values of this Lyapunov fleld identify regions 
of strong separation between particles, i.e., regions that will exhibit strong stretching 
during evolution, that can be identifled with the LCS [HI HOj [6]. 

In principle, to compute FSLE in 3d, the method of [6j can be extended to 
include the third dimension, by computing the time it takes for particles initially 
separated by do = [(Axq)^ + (A^o)^ + (Azq)^]"^^^ need to reach a flnal distance of 
df = [(Ax/)^ + (Ayj)^ + (Az/)^]-^/^. We will proceed this way for the turbulent channel, 
but, as indicated in vertical displacements are much smaller than horizontal ones 
in ocean flows. Therefore, the displacement in the z direction does not contribute 
signiflcatively to the calculation of rfj in the ocean, which prompt us to implement a 
quasi-3d computation of FSLEs: we use the full 3d velocity fleld for particle advection 
but particles are initialized in 2d horizontal ocean layers and the contribution Azf is 
not considered when computing df (see more details in [31j). In any case, since we allow 
the full 3d trajectories of particles, we take into account the vertical dynamics of the 
oceanic flows. 



Concerning the turbulent channel, where we can implement a fully 3d computation 
of the FSLE, we proceed as follows. A grid of initial locations xq = (xi^yj^Zk) is set 
up at time t, flxing the spatial resolution of the FSLE fleld (Fig. [T]). Particles are 
released from each grid point and their three-dimensional trajectories are calculated. 
The distances of each neighbor particle with respect to the central one (initially do) is 
monitored until one of the separations reaches a value df. 

In both systems considered, we obtain two diflFerent types of FSLE maps by 
integrating the three-dimensional particle trajectories backward and forward in time: 
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Figure 1. Computational setup for the calculation of the FSLE field in 3d. The 
FSLE at the location of the central particle (o) is a measure of the time it takes for any 
of the neighbor particles (•) to diverge from the central particle by a distance greater 
than 5f. 

the attracting LCSs (for the backward), and the repelling LCSs (forward) [IQIE]. We 
obtain in this way FSLE fields with a spatial resolution given by do. When a particle 
leaves the velocity field domain or reaches a no-slip boundary, the FSLE value at its 
initial position and initial time is set to zero. If the interparticle separation remains 
smaller than Sf past a maximum integration time At, then the FSLE for that location 
is also set to zero. 

2.2. Lagrangian Coherent Structures 

The identification of LCS calculated from Lyapunov fields in 2d fiows is straightforward 
since they practically coincide with (finite-time) stable and unstable manifolds of 
relevant hyperbolic structures in the fiow O |9l [10] (but see j32[ [IS]). The structure 
of these manifolds in 3d is generally much more complex than in 2d [23l |33] , and they 
can be locally either lines or surfaces. 

DiflFerently than 2d, where LCS can be visually identified as the maxima of the 
FSLE field, in 3d they are hidden within the volume data and one needs to explicitly 
compute and extract them, using the definition of LCSs as the ridges of the FSLE field. 
A ridge L is a co-dimension 1 orientable, differentiable manifold (which means that for 
a 3d domain ridges are surfaces) satisfying the following conditions [15j : 

(i) The field A attains a local extremum at L. 

(ii) The direction perpendicular to the ridge is the direction of fastest descent of A at 
L. 

This definition means that as we move away from the ridges, perpendicularly to it, we 
would see the fastest descent of the FSLE field. The method used to extract the ridges 
from the scalar field A(xo, t) is from [34j. It uses an earlier f35] definition of ridge in the 
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context of image analysis, as a generalized local maxima of scalar fields. For a scalar 
field / : ^ M with gradient g = V/ and Hessian H, a d-dimensional height ridge is 
given by the conditions 



where a^, i G {1, 2, . . . , n}, are the eigenvalues of H, ordered such that ai > . . . > a, 
and Bi is the eigenvector of H associated with ai. For n = 3, Eq. ^ becomes 



In other words, in the 61,62 eigenvectors point locally along the ridge and the 63 
eigenvector is orthogonal to it, so the the ridge maximizes the scalar field in the normal 
direction to it and in this direction the field is more convex than in any other direction, 
since the eigenvector associated with the most negative eigenvalue is oriented along the 
direction of maximum negative curvature of the scalar field. 

The extraction process progresses by calculating the points where the ridge 
conditions are verified and the ridge strength \as\ is higher than predefined threshold s 
so that ridge points whose value of is lower (in absolute value) than s are discarded 
from the extraction process. Since the ridges are constructed by triangulations of the set 
of extracted ridge points, the strength threshold greatly determines the size and shape 
of the extracted ridge, by filtering out regions of the ridge that have low strength. The 
reader is referred to [34] for details about the ridge extraction method. The height ridge 
definition has been used to extract LCS from FTLE fields in several works (see, among 
others, [36j). 

Since the A value of a point on the ridge and the ridge strength are only related 
through the expressions ^ and ([3]), the relationship between the two quantities is not 
direct, which makes difficult to choose the appropriate strength threshold s. A too 
small value of s will result in the extraction of very small LCSs that appear to have 
little inffuence on the dynamics, while a large value will result in only a partial rendering 
of the larger and more significant LCS, limiting the possibility of observing their real 
impact on the fiow. 

The ridges extracted from the backward FSLE map approximate the attracting 
LCSs, and the ridges extracted from the forward FSLE map approximate the repelling 
LCSs. The attracting ones are the more interesting from a physical point of view [6lfT2]. 
since particles (or any passive scalar driven by the fiow) typically approach them and 
spread along them, so that they are good candidates to be identified with the typical 
filamentary structures observed in tracer advection. 

3. Turbul6nt chann6l flow 

Turbulent channel fiow is a fiow between two stationary, parallel walls separated by a 
distance 25. It has been studied extensively due to its geometrical simplicity and its 
wall-bounded nature, which makes it a platform to study more complex turbulent shear 
fiows of greater technological interest. 



\/d < i < n , g^6^ = and ai < 0, 



(2) 



g^63 = and as < 0. 



(3) 
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The coordinates of the flow are: x for the streamwise direction, y for the cross- 
stream coordinate that separates the two plates, and z for the spanwise direction. The 
flow is maintained by a downstream pressure gradient ^ acting against the waU shear 
stress. The laminar flow solution Uq is a cross-stream parabolic proflle given by 

where /i is the dynamic viscosity. Following the Reynolds averaging method ^37j, the 
turbulent flow velocity u is decomposed in a mean U = ([/(^),0,0) and a fluctuating 
component = {u\v\w'). The mean turbulent velocity proflle U{y) diflFers from the 
laminar one, U{){y)^ by a lower centerline velocity t/(0) and increased near- wall velocity 
giving it a flatter shape. Due to the increase in mean velocity near the wall, the shear 
stress near the wall is higher for the turbulent case. The total shear stress r appearing 
in the averaged Reynolds equations gets contributions from both the viscous stress and 
the Reynolds stress —u'v' associated to the velocity fluctuations: 
T dU 



V— u'v' (5) 

P dy 

V = ji/ p the kinematic viscosity. The symmetries of the domain and the Reynolds 
equations imply that r depends only on the cross-stream coordinate ^, and the 
dependence is linear, so that it can be written as 

(6) 



The shear velocity Ur gives the velocity scale of the turbulent velocity fluctuations. The 
formula [37J : 



2 dU{y) 



(7) 

y=0 



dy 

allows to compute Ur from measurements of the mean velocity proflle from the 
simulations. A length scale can be formed by combining Ur with u: the wall scale 
6^ = u/ur- The wall distance can now be expressed as y^ = y/S^^ and the same 
normalization could be done for the rest of coordinates. The viscous Reynolds number 
Rcr = 5/5^ is simply the ratio between the two relevant length scales. 

The existence of coherent structures in turbulent wall-bounded flows has been 
known for several decades from investigations on intermittency in the interface between 
turbulent and potential flow regions, on the large eddy motions in the outer regions of 
the boundary layer, and on coherent features in the near- wall region ([38] and references 
therein). Since then, through experimental and numerical investigations, a picture of the 
organization of these coherent structures in the turbulent boundary layer has emerged, 
which has become rather complete from the Eulerian point of view [381 |39j. Our 
approach is a contribution to the Lagrangian exploration of these coherent structures, 
as in [22J. 

The longitudinal velocity fleld in the inner region of the channel (the viscous 
sublayer adjacent to the wall and the intermediate buffer region) is organized into 
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alternating streamwise streaks of high and low speed fluid. Turbulence production 
occurs mainly in the buffer region in association with intermittent and violent outward 
ejections of low-speed fluid and inrushes of high-speed fluid towards the wall. The outer 
region is characterized by the existence of three-dimensional (^-scale bulges that form on 
the turbulent/potential flows interface. Irrotational valleys appear at the edges of the 
bulges, entraining high-speed fluid into the turbulent inner region. A central element in 
the structure of the turbulent boundary layer is the hairpin vortex, mainly because it 
is a structure with the capability of transporting mass and momentum across the mean 
velocity gradient and because it provides a paradigm with which to explain several 
observations of wall turbulence [38l HQ] • 

3.1. Data 

The data used to extract the LCS come from a direct numerical simulation (DNS) 
of turbulent channel flow at a viscous Reynolds number Rcr = 180. The setup of 
the simulation follows that of [41j and is summarized in table [l] The simulations 
were conducted using the CFD solver Channelflow. The Channelf low code solves 
the incompressible Navier-Stokes equations in a rectangular box with dimensions 
Lx X 2S X with periodic boundary conditions in the x (so that fluid leaving the 
computational domain in the direction of the mean flow dit x = reenters it at x = 0) 
and in the spanwise z direction. No-slip conditions are imposed on ^ = ±6. The 
unsteady velocity fleld u is represented as a combination of Fourier modes in the x and 
z directions and of Chebyshev polynomials in the wall-normal direction. The pressure 
gradient necessary to balance the friction at the walls was chosen as to maintain a 
constant bulk velocity of |[/o. Time stepping is a 3rd-order Semi-implicit Backward 
Differentiation. Note that in our computations 6^ = l/Rcr = 0.0058 so that in wall 
units <y+ < 344. 

The flow was integrated from an initial base-flow with parabolic proflle and a 
small disturbance that evolved into a fully developed turbulent flow. The total 
integration time was At = 600 time units that in dimensionless form t^ = t{u^/v) 
gives At+ = 83.54. After an initial transient of about 200 time units the simulations 
reached a statistically stationary state from which statistics was accumulated. 

The mean quantities and flrst order statistics of our simulations where compared 
to those of j41j and the agreement is quite good. The proflle of the mean velocity in 
wall units is shown in Fig. [2j The proflle for the Reynolds stress —u'v' shows that the 
maximum (in absolute value) is located at approximately y+ = 30, in the outer limit of 
the buffer layer (sse Fig. [3]). 

3.2. Results 

The LCS were extracted from the turbulent velocity fleld data described in the previous 
section. A calculation of FSLE fleld in the entire turbulent channel was conducted in 
order to understand the statistical properties of the FSLE fleld in this class of turbulent 
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Figure 3. Reynolds stress u^v^ profile at Rcr = 180. Solid line: our simulations; 
squares: [41^ (given up to the channel centerline). 
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Table 1. Simulation parameters. Quantities with + refer to wall units. L^^ 25 and 
are the domain sizes in the y and z directions. Ax+, Ay~^ and Az~^ are the respective 
spatial resolutions (given at the first point above the wall for the y case), and n^^ riy 
and Uz the corresponding number of grid points. Re = US/u is the Reynolds number 
based on the channel center mean speed, whereas Rcr = u^d/v is the viscous Reynolds 
number. The nominal value is an input to the computer code, and the actual value 
comes by using Eq. ^ for the computed mean profile U{y). 



Re channel center 


4000 


Rcr nominal 


180 


Rcr actual 


172 


Lx 


47r 


S 


1 


Lz 


1- 


Lt 


2166.61 


(5+ 


0.0058 


Lt 


722.20 


rix 


128 


Uy 


129 


Uz 


128 




17.06 


Ay+ 


0.005 




5.6867 



Table 2. FSLE calculation parameters, dt is the integration time step and At the 
maximum integration time. 



Calculation 


do 


df/do 


At 


dt 


Complete channel 


0.024 


20 


172 


0.05 


LCS subdomain 


0.003 


67 


10 


0.05 



flows. A subsequent calculation in a subdomain of the channel was used to extract the 
LCS in that subdomain for a sequence of time instants. The setup of both calculations 
is shown in table O 

3.2.1. The 3d FSLE field. The 3d backward FSLE field for the entire channel was 
calculated at a single time instant in the statistically steady state. The initial and final 
distances do and df were chosen as a balance between encompassing the widest possible 
range of scales of motion (measured by the ratio rf//rfo), and adequate resolution and 
computational cost. The initial distance is of the order of 4:6^ and the final distance of 
the order of 0.56 so that the ratio of scales, df/do^ is approximately i?e^/8. 

Figure [4] shows an instantaneous configuration of the FSLE values in a 
streamwise/ wall-normal plane. Maxima of the FSLE field appear to be organized 
into sloping structures located in the region 20 < < 100. This organization bears 
resemblance to the widely accepted picture of organized structures in wall turbulence 
where the outer region is dominated by packets of sloping hairpin vortices [iQl |38] . The 
channel center is devoid of high FSLE values but coherent patches of low FSLE values 
can still be observed. 

A cross-stream FSLE profile is obtained by averaging the 3d field over the periodic 
directions x and z. It is shown in Fig. [5| The profile is symmetric about the channel 
centerline and shows a maximum at approximately y+ = 30, at the same location 
of the maximum in the Reynolds stress —u'v^. Because of the periodic boundary 
conditions in the x and z directions the average profiles along these directions are 
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500 1000 1500 2000 

x+ 

Figure 4. FSLE values shown on a streamwise/wall-normal plane in the turbulent 
channel. Walls are at the top and bottom of the figure. Mean velocity is towards the 
right. 



rather unstructured, and we resort to two-point correlation functions to quantify the 
statistical structure properties. For each plane parallel to the walls, i.e. for each value 
of we compute the fluctuations of the FSLE values around the average in that plane: 
A(x+,^+,z+) = A(x+,^+,z+) — (A(x+, From this quantity we define the 

streamwise Rxx{y^]x^) correlation function as: 

^ , + (A(x+,2/+,z+)A(x+ + x+,2/+,z+))^+^^+ 

Rxx[y^] x^) = /w , , ,X2\ — 8 

(A(x+,2/+,z+)2)^+^^+ 
and the spanwise Rzz{y^] z^) correlation function 

^ , + (A(x+,2/+,z+)A(x+,?/+,z+ + z+))^+^+ 

(A(x+,2/+,z+)2)^+^^+ 

In the above equations the averages are over the periodic directions x+ and z^. The 
correlations are shown in Figs. [6] and [7] at different distances from the walls: one 
smaller, one larger, and one approximately coincident with the location of the maximum 
Reynolds stress. These functions reveal sizes and organization of the different structures 
in the Lagrangian FSLE field, to be contrasted with Eulerian correlation functions in 
the same system 
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Figure 5. FSLE profile averaged over as a function of the cross-stream 

normalized coordinate . Only half of the channel is shown since the profile is quasi- 
symmetric about the channel center line. 

Close to the wall {y^ = 12.2), viscous effects make the structures smooth and 
extended along the streamwise direction (see the long correlation length for this case 
in Fig. [g]), forming streaks. In the transverse direction the oscillations seen in R^z 
for = 12.2 indicate an approximately periodic arrangement of the streaks [22], with 
a spacing in the range 100 — 150 wall units. This pattern of organization is similar to 
what is seen in Eulerian descriptions 021 [38]. At planes further away from the wall 
(y^ = 28.4 and y^ = 122.1 in Figs. |6]and[7|, correlation functions in both directions 
become shorter ranged, and periodic features are progressively lost. This corresponds 
to a rather disordered distribution of structures, each with a typical size related to the 
width of the correlation functions, i.e. of the order of 50 wall units, as also seen in Fig. 

SI 

3.2.2. The 3d LCS. The previous description summarized the statistical properties 
of the different structures appearing in an instantaneous FSLE field. To make further 
progress we now extract three-dimensional attracting LCSs in a region of the channel 
at a series of time instants. The extraction domain had dimensions Lt x Lt x Lt = 
103 X 129 X 124. The initial separation do and distance ratio df/do were increased from 
the previous calculation to improve the resolution and extract smoother structures, but 
sacrificing a complete view of 3d LCS in the turbulent channel. The extraction threshold 
was set to s = 50000, a compromise value between speed and cost of extraction and 
continuity of the extracted surfaces. The FSLE fields were calculated for an interval of 
1.5 time units with a time step of 0.1 units. 
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_0.2i ^ ^ ^ ^ ^ ^ ' 

50 100 150 200 250 300 350 

Figure 7. Spanwise correlation function Rzz{y^]z'^) as a function of the spanwise 
separation at three distances from the lower wall: Continuous black line: = 12.2; 
blue dashed line = 28.4; red dot-dashed line: y^ = 122.1. 
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Figure 8. 3d attracting LCS in the channel flow together with a FSLE map at the 
fixed plane x = 6.0. Time goes from top to bottom, at intervals of 0.1 time units. The 
fiow direction is in the positive x direction in each panel, and a wall is at the bottom. 
The sequence shows how one of the flow structures is advected and passes through the 
X = 6.0 plane. 
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3d LCSs are rendered in Fig. [8| in a sequence of time instants, as they pass through 
the calculation domain. They have a clearly 3d shape and move with the flow. The LCS 
seem to create a boundary between the inner turbulent region and the outer region that 
is practically devoid of FSLE. The highest LCS have (^-scale heights above the wall, and 
have a distinct mushroom shape enclosing the regions of the channel closer to the wall, 
where high FSLE values can be found. Near the wall, the LCS adopt the shape of sheets 
parallel to it, which reflects the high rates of shear that occur in that region. These 
sheets form the base of the mushroom-shaped excursions up to the channel center. 

4. Oceanic flow 

Contrarily to the turbulent flow of the previous section, large scale oceanic flows, 
naturally turbulent, can be considered as almost 2d due to rotation and stratiflcation 
effects. This fact makes the theory of 2d turbulence a very important tool to understand 
the ocean processes that occur at large scales. The main characteristic of 2d turbulence 
is the existence of an inverse energy cascade, from the small to the large scales and a 
direct enstrophy cascade. This cascade manifests itself by the creation of large coherent 
vortices, and by the process of fllamentation by which strain regions in the boundaries 
of the vortices produce lines of vorticity that are continuously stretched and deformed 
by the flow, concentrating the vorticity gradient in the small scales. This behavior is 
often observed in oceanic flows thereby conflrming the importance of the 2d turbulent 
processes. 

4.1. Data 

The Benguela ocean region (Fig. [9]) is situated off the west coast of southern Africa. It 
is characterized by a substantial mesoscale activity in the form of eddies and fllaments, 
and also by the northward drift of Agulhas eddies. 

The velocity data set comes from a regional ocean model (ROMS) simulation of 
the Benguela Region [43j. ROMS [SJ US] is a split-explicit, free-surface, topography 
following model. It solves the incompressible primitive equations using the Boussinesq 
and hydrostatic approximations. Potential temperature and salinity transport are 
included by coupling advection/diffusion schemes for these variables. The model was 
forced with climatological data. The data set area extends from 12°S to 35°S and from 
4°E to 19°E (see Fig. |9]). The velocity fleld u = (u^v^w) consists of two years of daily 
averaged zonal (u)^ meridional (v)^ and vertical velocity (w) components, stored in a 
three-dimensional grid with an horizontal resolution of 1/12 degrees (^ 8 km), and 32 
vertical terrain-following levels. Additional details can be found in [3T] . 
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Figure 9. Benguela ocean region. Grey region is Southwestern Africa. The velocity 
field domain is limited by the continuous black line. The FSLE calculation area is 
limited by the dash-dot black line. Bathymetric contour lines are from ETOPOl 
global relief model [46] from m depth up to 4000 m at intervals of 500 m. 



4^2. Results 

4.3. Three-dimensional FSLE field 

The three-dimensional FSLE fields were calculated for a 30 day period starting 
September 17 of year 9, with snapshots taken every 2 days. The fields were calculated 
for an area of the Benguela ocean region between latitudes 20°S and 30°S and longitudes 
8°E to 16°E (see Fig. |9]). The calculation domain extended vertically from 20 up to 580 
m of depth. Both backward and forward calculations were made in order to extract the 
attracting and repelling LCS. 

Figure [To] displays the vertical profile of the average FSLE for the 30 day period. 
The small differences between the backward and the forward values are due to the 
different intervals of time involved in their calculation. There is a general decrease with 
depth, with a notable peak in the profiles at about 100 m. The reason for this local 
maximum in the FSLE profiles is not clear but it could be due to increased vertical 
shear enhancing the mixing rates at those depths [31j . 



In the left panel of Fig. [TT] a snapshot of the attracting LCSs for day 1 of the 
calculation period is shown. The structures appear as thin vertical curtains, most of 
them extending throughout the whole depth of the calculation domain. The horizontal 
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Figure 10. Vertical profile of a 30 day average of backward (squares) and forward 
(circles) FSLE. The 30 day average 3d field was further spatially averaged over each 
horizontal layer to produce the vertical profiles. 



slices of the FSLE field in Fig. [IT] (left panel) show that the attracting LCS fall on 



the maximum FSLE field lines. The atracting and repelling LCS (Fig. 11, right panel) 
populate the calculation region, testifying the enhanced mixing activity that is known to 
occur in that particular ocean region. The quite entangled "web" in which attracting and 
repelling LCSs intersect mutually provides the skeleton for the barriers and pathways 
controlling transport [6l [H] . 

At this point, it may help to stress the differences between the Eulerian and 



Lagrangian detection of coherent structures. This can be seen in Fig. [12] where the 
boundaries of a mesoscale eddy are shown using the Q-criterion and the attracting and 
repelling LCS. The Q-criterion [47] uses the second invariant of Vu: 



Q = li 



lOI 



|S||^), 



(10) 



where = tr(nri^), ||S|p = tr(SS^), and 12, S are the antisymmetric and 

symmetric components of Vu, to identify regions where rotation dominates strain 
(Q > 0), commonly identified with coherent vortices, and strain dominated regions 
{Q < 0). We refer the reader to [48j and [49j for reviews and criticism of several 
Eulerian criteria. 

Eulerian and Lagrangian measures limit approximately the same region, but are 
substantially different. The Q-criterion is related to the instantaneous configuration of 
the second invariant of Vu and therefore conveys only local information about fluid flow 
processes. The Lagrangian perspective, on the other hand, provides an integration of 
the temporal evolution of material properties of the flow, e.g. material transport, and 
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Figure 11. 3d LCS in the Benguela region for day 1 of the calculation period. Left 
panel: Attracting LCS together with horizontal slices of the backward FSLE field at 
120 m and 300 m depth. Right panel: Attracting (blue) and repelling (green) LCS. 
Colorbar refers to colormap of horizontal slices in the left panel. The units of the 
colorbar are day~^. 



thus should give more meaningful information about the processes that rely on ocean 
material transport. 

This issue can be further explored by looking at a filamentation event (described 
more extensively in [S]). A set of particles were released inside the eddy at day 1 at 



a depth of 50 m. At day 11 of the calculation period (see Fig. 13), they have formed 



a filament that is expelled from the eddy, so that particles clearly cross the Q-criterion 
isosurface. This shows that the Eulerian criteria is inadequate as an indicator of regions 
of material transport in the flow. On the contrary, it can be observed that the Lagrangian 
description of the eddy boundaries does bear relation with material transport into and 
out of the eddy, since the particle filament leaves the enclosed region that we associate 
with the mesoscale eddy by following one of the identified Lagrangian boundaries. 



5. Conclusions 



Lyapunov exponents are useful to identify Lagrangian Coherent Structures in turbulent 
fiows. These constitute the pattern determining the pathways of particle transport in 
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Figure 12. Attracing (blue) and repelling (green) LCS on day 1 of the calculation 
period together with Q-criterion isosurface at Q = 10~^^ (red). 

the flow. They strongly influence the transport and mixing properties in the fluid. 

In this paper we have used a particular type of Lyapunov exponents, the so-called 
Finite-Size Lyapunov exponents, to identify LCS in 3d flows. The flnite size Lyapunov 
exponent was used to measure the rate of streching of initially nearby fluid particles 
in the flow domain and the Lagrangian coherent structures where identifled as the the 
ridges of the FSLE fleld. These ridges were flltered in order to retain only the strongest 
attracting or repelling structures. 

In a turbulent channel flow, the LCSs appear as mushroom-shaped excursions of 
near- wall sheet-like structures of a scale comparable to the channel width. They separate 
the channel into an interior region, where the FSLE attains high values, and an exterior 
region, showing low FSLE values. The distribution of LCS in the turbulent channel 
resembles the commonly accepted picture where upward excursions of near wall fluid 
coexist with inward rushes of mid-channel irrotational flow. Further work is necessary 
to elucidate the relations between LCS and fluid transport in these type of flows, not 
least because the visualization of 3d structures and transport in turbulence is a complex 
and time-consuming subject. 

In a quasi-2d mesoscale oceanic flow, the LCSs appear as quasi-vertical surfaces 
highlighting the fact that dispersion in this case is mainly horizontal. The high mixing 
activity can be deduced from the proliferation of LCS in the flow domain and their 
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Figure 13. Attracing (blue) and repelling (green) LCS on day 11 of the calculation 
period together with Q-criterion isosurface at Q = 10~^^ (red). The particles (black 
dots) were released inside the eddy at day 1 at a depth of 50 m and are leaving now 
the eddy as a filament along the upper part of the attracting LCS. 



mutual intersection. These LCS were seen to provide barriers and pathways to transport 
in the case of a mesoscale eddy, contrary to Eulerian measures that failed to provide 
indicative locations or directions of major transport events. 

The results shown in this paper highlight the usefulness of Lyapunov analysis and 
dynamical systems theory as a tool to study transport and mixing in fluid flows, through 
the concept of Lagrangian coherent structures. 
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